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Abstract 

Randomly packing spheres of equal size into a container consistently results in a 
static configuration with a density of ~64%. The ubiquity of random close packing 
(RCP) rather than the optimal crystalline array at 74% begs the question of the 
physical law behind this empirically deduced state. Indeed, there is no signature of any 
macroscopic quantity with a discontinuity associated with the observed packing limit. 
Here we show that RCP can be interpreted as a manifestation of a thermodynamic 
singularity, which defines it as the "freezing point" in a first-order phase transition 
between ordered and disordered packing phases. Despite the athermal nature of 
granular matter, we show the thermodynamic character of the transition in that it is 
accompanied by sharp discontinuities in volume and entropy. This occurs at a critical 
compactivity, which is the intensive variable that plays the role of temperature in 
granular matter. Our results predict the experimental conditions necessary for the 
formation of a jammed crystal by calculating an analogue of the "entropy of fusion" . 
This approach is useful since it maps out-of-equilibrium problems in complex systems 
onto simpler established frameworks in statistical mechanics. 
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Since the time of Kepler it is thought that the most efficient packing of monodisperse 
spherical grains is the face centered cubic (FCC) arrangement with a density of 74 % 
Thus, we might expect that spherical particles will tend to optimize the space they occupy 
by crystallizing up to this limiting density. Instead, granular systems of spheres arrest in a 
random close packing (RCP), which is not optimal but occupies ~64% of space Previous 
studies have derived geometric statistical models to map the microscopic origin of the much 

nn 

debated 64% density of RCP [2—9] . However, the physical laws that govern its creation and 
render it the most favorable state for randomly packed particles remains one of the most 
salient questions in understanding all of jammed matter js-5, 8]. For instance, while it is 
known that systems in equilibrium follow energy minimization and entropy maximization 
to reach a steady state, the mechanism by which RCP is achieved is much sought after. 

Here we propose a thermodynamic view of the sphere packing problem where the exper- 
imentally observed RCP can be viewed as a manifestation of a singularity in a first-order 
phase transition. Despite the inherent out-of-equilibrium nature of granular matter, the 
formation of a jammed crystal can be mapped to a thermodynamic process that occurs at 
a precise compactivity where the volume and entropy are discontinuous. 

We investigate mechanically stable packings ranging from the lowest possible volume 
fraction of random loose packing (RLP) n]] to FCC. We numerically generate packings of 
N =10,000 spherical particles of radius R = lOOum in a periodically repeated cube. Initially, 
we use the Lubachevsky-Stillinger (LS) 111, [12] and force-biased (FBA) 13] algorithms to 
generate amorphous configurations of unjammed hard-spheres fluids at infinite kinetic pres- 
sure and volume fraction (see Appendix-Section fl}. While these configurations are 
geometrically jammed, they are not jammed in a mechanical sense since the particles do not 
carry any forces: the confining stress a (not kinetic) is zero. In fact a key difference between 
granular materials jammed under external stress or gravity and hard-sphere fluids is that, in 
the former, each particle satisfies force and torque balance. In order to study mechanically 
stable packings characterized by a jamming transition we introduce interparticle forces via 
the Hertz-Mindlin model of normal and tangential forces allowing the particles to be soft 
but with a large Young modulus, Y. We then simulate the process of jamming by Molecular 
Dynamics (MD) simulations using previously developed methods 14] to compress the LS 
and FBA packings from <pi to a final jamming density, <pj. Ultimately, we obtain mechan- 
ically stable packings just above the jamming transition (in the limit of vanishingly small 
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confining stress, or equivalently in the hard-sphere limit, a/Y — > + ) covering a range of 4>j 
from r ip = 0.55 to crystallization at 0f cc = 0.74. 

The mechanical coordination number averaged over all the particles in a packing, Zj, 
characterizes different states of granular matter 0, Q, [l5]. Therefore, our study begins 
by plotting Zj versus 4>j f° r & U generated packings. Figure [T] suggests the existence of a 
transition occurring at RCP evidenced by the abrupt plateau in Zj. This transition could 



be thought of as an ana 
thermal equilibrium 16 



ogue to the classical hard sphere liquid-solid phase transition in 



17| . Such an analogy becomes apparent if one identifies Zj of the 



jammed packing with the kinematic pressure of the equilibrium hard sphere system [6[ and 
it is in agreement with a recent conjecture regarding the definition of RCP js]. 

Figure [T] identifies two branches and a coexistence region: (i) An ordered branch of 
crystallized states with <pj ranging from 0.68 to a FCC lattice at 0.74. (ii) A disordered 
branch within 0.55 ~ 0.64 which can be fitted with the statistical theory of [7j]: 4>j = Zj/ (Zj + 
2y / 3) as shown in the figure. (Hi) A coexistence region between 0.64 to 0.68 displaying a 
plateau at the isostatic coordination number, Zi SO = 6 j^, 14, 18]. The intersection between 
the disordered branch and the coexistence line identifies the "freezing point" of the transition 
providing a definition of RCP. Using the theoretical results of |7j, freezing occurs at Z iso = 6 
and rcp = 6/(6 + 2\/3) ~ 0.634. The corresponding "melting point" appears at the other 
end of the coexistence at me it = 0.68, signaling the beginning of the ordered branch. Finite 
size analysis is shown in the Appendix-Fig. [5]A the results for 500 and 10,000 spheres 
are consistent with each other. Other geometric aspects of the transition are discussed in 
Appendix-Section [III 

To reveal the nature of the newly found phases we start with a descriptive viewpoint and 
then turn to a thermodynamic analysis to model the transition. In order to investigate if 
the concept of phase transition applies to the trend observed in Zj, one commonly looks at 
the global (Qi,Wi), and local (qg) orientational order parameters for a signature of varying 
amounts of crystallization present in the packings as defined elsewhere 19( (see Appendix- 
Section [TTlA] and Fig. [2] for definitions). The salient feature of Qi is that its zero value 
means disorder and non-zero value means crystallization. Therefore, the presence in Fig. 
[2A of an increase in Qi from zero at rcp defines the beginning of the coexistence reg ion. 
Typically, a first-order transition is marked by a nonzero third-order invariant W\ |l ol . I20I ] 
which we find appears at the melting point me it signaling the onset of the ordered branch 
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(Fig. W). 

Interestingly, the ordered phase has two significant peaks in the probability distribution 
P(g 6 ) of the local order parameter of each particle, q 6 (Fig. |2C). The peaks correspond to 
FCC and HCP [19] signaling that both crystalline configurations are present in the ordered 
structure. From the available data we cannot rule out the possibility of another transition 
from HCP to FCC before <pj ~ 0.74. The Gaussian distributions P{qe) obtained for <pj within 
0.55 ~ 0.64 show no preferred lattice structure in the random branch. While the relative 
peak positions in P(q$) do not change, the percentage of crystal and random phase found in 
the packing progresses from one to the other in the coexistence region. Microscopically, the 
existence of the two pure phases is starkly represented by the two separated distributions of 
local Voronoi volume fractions P(4> vor ) for which the same phenomenology of P(qe) applies 
(Fig. [2p, (j) VOI = Vg/V V oT, where V vor is the Voronoi volume |7|, [l5| of each particle of volume 
V g ). This descriptive analysis is further supported in Appendix-Section IIHI 

Having identified the structure of the phases we now develop a thermodynamic viewpoint 
of the RCP transition to rationalize the obtained results. Transitions in equilibrium physical 
systems are driven by a competition of energy and entropy. Instead, a transition in athermal 
jammed matter is driven by the minimization of the system's volume W by compactification 



and entropy maximization of jammed configurations [2ll. |22|. In accordance with the second 
law of thermodynamics, the granular system tends to minimize its Gibbs-Helmholtz "free 
energy" F = W — XS rather than W alone. The compactivity of the system X = dW/dS 
is a measure of how much further compaction a packing can undergo; the lower the volume 
the lower the compactivity 21] . Thus, we map the packing problem to a thermodynamic 
problem where the volume W replaces the energy and X takes the role of temperature. The 
principle of free energy minimization can thus be applied. 

The equations of state, <fij{X) and S(<f>j), can be calculated from the fluctuations of the 
Voronoi volumes in the disordered and ordered phases (23 1 , cri(<f>j) and ^(^j) respectively, 
in analogy to the standard Boltzmann statistical mechanics (a 2 = (w^ OT ) — (w VOI ) 2 and 
Wvor = l/0vor = V VOT /V g is the reduced Voronoi volume). Figure EJ\ shows clearly the 
existence of the two pure phases and a discontinuity between both branches. We obtain 



the compactivity by integration of <j\ and a 2 using Einstein fluctuation theory 23N25| (see 



Appendix-Section IIV Al for more details, we set ks — 1 for simplicity, X is given in units of 
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Vg and entropy is dimensionless) : 



dtp 1 

1 1 [+' d<p 1 , , 

where A" me it = X(0 me i t ) and X r i p = AT(0 r ip) are the compactivities at the melting point and 
at RLP, respectively. Once XUj)j) is obtained from Eq. ((T]), the entropy density, s = S/N, 
is calculated by integration [23| (see Appendix-Section IIV Aj) : 



'j) = s rcp + V g — , rip < 0j < rcp , (2a) 



A 

»(<>.,) = }'., I x(^--- ' ° mAl ~ °' ~ (21>) 



where we have used that the entropy of FCC is zero in the thermodynamic limit. There 
are three unknown integration constants in Eqs. and (jSJ): A" r i p , X mc i t and the entropy 
of RCP, s rcp . To close the system, we consider the conditions for equilibrium between 
the phases [21]: (a) "thermal" equilibrium X me i t = X vcp = X c , where X c is the critical 
compactivity at the transition, and (b) the equality of the free energy density (or chemical 
potential), f = F/N, at the melting and the freezing RCP point: / me i t = f Tcp . This implies, 
uj Tcp -(X c /Vg)s rcp = u melt -(X c /V g ) s me it, where w = l/cpj = W/(NV g ) is the reduced volume 
of the system. The third integration constant X rlp can be considered infinite [7|, |23j, [25] 
since RLP is the highest volume of the system. While the precise value of X T \ P does not 
affect our conclusions, a more accurate finite value can be obtained by fitting the entropy 
Eq. d2J) with an independent measure obtained by cluster analysis from information theory 
(Shannon entropy, s sria n) as developed in [23] (Appendix-Section IIV Bp . Figure [3]B shows 
that the entropy from the thermodynamic integration Eq. (|2j) and s s han agree well (up to a 
multiplicative constant) supporting the framework of Eqs. ([I])-(j2j). The entropy is composed 
of two branches plus the coexistence region (green line in Fig. 03). 

FigureHJA. displays a discontinuity in s(X) at X c = 0.031K, revealing the first-order nature 
of the transition which is accompanied by an "entropy of fusion" Asf us = s rcp — s me \ t = 3.0. 
The volume fraction is discontinuous at X c (Fig. 33) where the system jumps from RCP to 
the melting point releasing an amount of volume given by the "enthalpy of fusion" A/if us = 
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X c Asf us = 0.09V g while the compactivity stays constant. This process corresponds to the 
typical latent heat in exothermic first-order transitions. 

Systems jammed at RCP need to overcome a volume barrier Ahf ns for crystal formation 
or, equivalently, particle displacements Arf us w 0.45-R. From a thermodynamic perspective, 
the requirement is equivalent to bringing a random packing in contact with a compactivity 
bath at X < X c = 0.031V^. The fundamental idea is to surround a random packing 
above X c with a crystal lattice below X c and perturb the system to equilibrate. A shear 
cycling experiment — which conserves the shape of the box containing the particles — suffices 
;o explore the crystal branch (3, 0|- Shear-induced crystallization has been observed 



26l . 1271 ] when the maximum angle of horizontal shear is above 9 ~ 10°. This value is of 



the same order as our estimate of the shear amplitude to crystallize at <pj = 68% based on 
the entropy of fusion, which gives 9 = tan _1 (Arf us /2i?) « 13°. Furthermore, recent shear 
cycling experiments js), S| appear to be in reasonable agreement with the present results. 
We also expect that 2d equal-sized disks may have a near zero entropy of fusion owning to 
their tendency to easily crystallize while Asf us may sharply increase in 4d and above [6]. 

The behavior of the free energy density shown in Figs. HP and HP summarizes the 
mechanism to achieve RCP. The free energy in Fig. HP increases as X decreases from RLP 
to freezing at RCP. At X c , the system transitions to the phase with the lower free energy 
through an entropy discontinuity given by s = —df/dX. The system may also enter the 
metastable branch as indicated in Fig. HP and in Figs. [T]and|3B from a — > b. In the spirit of 
Landau mean field theory of phase transitions, we relate the distribution of the local order 
parameter to the free energy functional, J 7 , and X as P(qe) ~ exp[— .F(g 6 )/Jf] [20]. Figure 
HP shows J-'i^qe) displaying the minima of J 7 ^) defining the order and disorder phases at 
different <ftj. We find that the location of the minimum at g™ in « 0.425 remains constant 
from RLP up to the freezing point as expected in the disordered phase. The value of J 7 ^™ 111 ) 
is very deep for = 0.55 and becomes less deep as the freezing point is approached. The 
value of J-"(g™ m ) at the freezing and melting points become approximately equally deep 
indicating the phase coexistence at X c . Within a statistical mechanics framework, these 
results are a natural consequence and give support to such an underlying statistical picture. 

In conclusion, treating granular packings from the perspective of theoretical physics de- 
veloped by Boltzmann and Gibbs has the potential to answer basic questions in the field of 
disordered media. State variables like the compactivity can be introduced with the potential 
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of identifying transition points between different phases, a fact that can unequivocally define 
RCP as the freezing point in a discontinuous transition. This formalism may be useful in 
analyzing other related_transitions in complex systems ranging from optimization problems 



m computer science 



29] to the physics of glasses [6]. Other unsolved packing problems in- 



cluding finding the densest arrangement of rods, ellipsoids, spherocylinders, binary mixtures, 
Platonic and Archimedean solids — which are known to pack better than spheres |9j, |30] — can 
now be analyzed from the proposed thermodynamic view of phase transitions. 
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FIG. [TJ The RCP transition. We plot the mechanical coordination number Zj versus 
the volume fraction <ftj for each packing. We identify: (i) a disordered branch which can 
be fitted with the statistical model of [t| as shown, (ii) a coexistence region, and (Hi) 
an ordered branch. Error bars are calculated over 523 packings obtained from initial LS 
configurations. The 3d plots visualize how the transition occurs in terms of arrangements 
of contacting particles. White particles are random clusters, light blue are HCP and green 
are FCC clusters. Further microscopic information regarding the transition is provided in 
Appendix-Section III Bt 

FIG. [2J Descriptive viewpoint of the RCP transition. (A) Global orientational order 
parameters Qi versus <pj for different packings signaling the freezing point at rcp . A linear 
fit is possible in the coexistence region. (B) Global third-order invariants W\ versus <ftj 
signaling the melting point at <p me \t- (C) Probability distribution of local orientational 
order parameter P{q%) versus q§ (vertical axis) for packings with (fcj (horizontal axis). For 
packings with <pj within 0.68 ~ 0.72, the distributions have two significant peaks centered at 



4« 



fee _ n K7 n„A „ hc P 



0.57 and q 6 = 0.48, which correspond to FCC and HCP arrangements, respectively 
19j . Color bar indicates the values of P{q§)- Since the peaks are very pronounced, we 
plot P{qo) up to the indicated value. (D) Probability distribution of local volume fractions 
of the Voronoi volumes of each particle, P(4> VO r) versus vor (vertical axis) for different <pj 
(horizontal axis). The plot indicates a clear discontinuity between both branches. Color bar 
indicates the values of P((f) VO r) which are plotted up to the indicated value. 

FIG. [31 Equations of state of the RCP transition. (A) Volume fluctuations of the Voronoi 
cell of a particle as a function of <pj. The data indicates a discontinuity between the ordered 
and disordered branches which are fitted by functions as indicated. These fittings are used 
in the integrations of Eq. ([1]). The larger fluctuations in volume observed in the order 
state compared to the disorder state at similar are due to the fact that the system packs 
better in the former and thus displays larger fluctuations when the system volume is the 
same. (B) Entropy obtained from fluctuation theory in Eq. ([2]), s, and Shannon entropy 
from information theory, s s han, versus volume fraction <pj. Both entropies agree (up to a 
multiplicative constant, k = 0.1, as indicated) confirming our calculations. The extended 
branch denotes a metastable state ending at point b at an hypothetical Kauzmann density, 
<Pk, in analogy with the physics of glasses [6J (see Appendix-Section III CI) . 

FIG. HI Thermodynamic viewpoint of the RCP transition. All the observables are con- 
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sistent with a transition at X c = 0.031V^. (A) Entropy s versus X. (B) Volume fraction 
4>j versus X. (C) Free energy density / versus X. We extend / for both branches to indi- 
cate the possible metastable states. At X c the system follows the minimization of the free 
energy signaling the transition from RCP to order. (D) Free energy functional J^iqe) versus 
q 6 (vertical axis) (fij (horizontal axis). Color bar indicates the values of ^(qe), which are 
plotted in the range indicated to focus on the region of coexistence. The minima correspond 
to the disordered phase and the FCC and HCP phases in the ordered branch. 
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Appendix 

A first order phase transition at the random close packing of hard spheres 

Yuliang Jin and Hernan A. Makse 



Here, we describe the details of the MD simulations (Section HJ), geometrical interpreta- 
tions of the transition (Section [TXD , and the calculations leading to the descriptive (Section 
TTT|) and the thermodynamic (Section IIVI) view of the RCP transition. 



I. ALGORITHM 

We use computer simulations to obtain jammed packings containing N = 10, 000 mono- 
disperse spheres of radius R = 100/im with periodic boundary conditions. We first apply 
a modified Lubachevsky-Stillinger (LS) algorithm ll|, Q to generate packings of densities 
up to ~0.72. In the algorithm, a set of random distributed points grow into nonoverlaping 
spheres at a fixed expansion rate 7. The spheres are considered perfectly elastic and evolve 
in time according to Newtonian dynamics. The configurations eventually arrive at out-of- 
equilibrium states with a diverging collision rate and a density 0j. Practically, we set the 
reduced kinetic pressure of the fluid defined as p = PV/NksT to be 10 12 [12|] as a criteria of 
the diverging collision rate. The final packing configurations depend on the expansion rate 
7: large values of expansion rate result in random packings with very low packing densities, 
while small values of expansion rate result in packings with higher densities. 

Although the packings obtained from the modified LS algorithm are considered as "ge- 
ometrically jammed", they are not jammed in the mechanical sense since the particles do 
not carry on any forces. In order to study mechanical stable packings characterized by a 
jamming transition, we model the microscopic interaction between deformable grains by the 



nonlinear Hertz-Mindlin normal and tangential forces [?], 



JjJ]. We use configurations from 



the modified LS as the starting point, (pi, and apply molecular dynamics to simulate Newton 
equations for the evolution of the particles following algorithms in 

an 

The aim of this part of the protocol is to generate mechanically stable jammed packings 
at the jamming transition <f)j. For the packings obtained by the LS algorithm, we first reset 
the velocities of the particles to zero. At this point the system stress a and mechanical 
coordination number Zj are zero since there is no deformation or overlapping between the 

16 



particles. We notice that the stress a is not the kinetic pressure, p, measured in the LS 
packings which diverges at the end of the LS protocol. Here a refers to the mechanical 
pressure related to the trace of the stress tensor via a = <Tu/3, where 



R 




(3) 




contacts 



where the sum is over all the contact forces, fl denotes the i-th component of the contact 
force, n c is the unit vector joining the center of two spheres of radius R in contact and V is 
the system volume. 

The system is then compressed isotropically by a constant compression rate until a given 
nonzero stress a is reached. Next, we turn off the compression and allow the system to relax 
with constant volume. If the system eventually reaches a jammed state with a fixed nonzero 
o and coordination number, the system pressure will remain unchanged over a large period 
of time (usually ~ 10 7 MD steps); otherwise, if the system is not stable, the pressure will 
relax to zero very fast Q. 

Previous studies show that there exists a jamming transition for granular matter as 



Here, <jj is zero for frictionless packings. However, it could has nonzero value for frictional 
packings, in general. The exponent a = 3/2 is trivially related to the Hertz-law of in- 
terparticle contact force and /3 = 1/2 seems to be universal over different force laws ^]. 

In practice, it is difficult to reach a jammed packing exactly at the transition point <ftj 
while it is much easier to get a stable packing with slightly higher pressure. In order to 
approach the transition point, a jammed packing at higher pressure than Oj in Eq. (jlj) 
obtained using the above protocol is decompressed with a negative compression rate until 
certain lower pressure is reached. Then the system is allowed to relax again to check for 
mechanical stability. If it is stable, then the system is decompressed further to an even lower 
pressure, and we check its stability again. By this process (called the split algorithm in {?|) 
we are able to approach the density <pj at the jamming transition point as close as possible 




(4) 



and 




(5) 
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FIG. 5: (A) The mechanical coordination number Zj versus the volume fraction <pj with different 
system sizes (500 and 10,000 particles) and algorithms (Lubachevsky-Stillinger LS, force-biased 
FBA, and split SA algorithms). The results show that the transition does not depend on the 
system size and the algorithm used. (B) Comparison between geometrical coordination number 
z g and mechanical coordination number Zj. Along the disordered branch, z g closely follows the 
mechanical coordination number Zj. We expect that in the thermodynamical limit the gap between 
both coordinations may diminish. z g and Zj start to diverge at RCP, which is an indication of 
increasing geometric degeneracies in the contact network at the onset of crystallization. 

within the system error. The pressure of the packings at the jamming point studied in this 
paper is 100 ± 8 KPa. The difference between the volume fraction of these packings and 
the critical volume fraction <pj from power-law fitting in Eq. fll]) is about 10~ 3 . 

The same preparation protocol is repeated by using the force-biased algorithm (FBA) of 



131 ] as initial protocol. The force-biased algorithm is a variant of the original method of 



W. S. Jodrey and E. M. Tory, Phys. Rev. A 32, 2347 (1985), We also generate packings 



ft 



0.3 



following the split algorithm of [7J starting with low initial volume fractions at fa 
below RLP. 

The mechanical coordination number, Zj, versus the final jamming density, fa, is plotted 
in Fig. [5J4. for all the obtained jammed packings which total 720. The plot signals the 



existence of a transition at RCP. It is analogous to the equilibrium 
hard-spheres, if we replace Zj by the kinetic pressure of the fluid |6|. 



iquid-solid transition in 
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II. OTHER ASPECTS OF THE RCP TRANSITION 



A. Finite size analysis 

It is important to determine the finite size effects of our results. Figure 0\ shows the 
results for smaller systems of 500 particles compared with 10,000 spheres system used in 
Fig. [TJ We find that both plots are consistent with each other. While Fig. [1] shows the 
average of Zj over the LS packings, Fig. shows each point representing a single packing 
obtained with the indicated algorithms for different system sizes. We find that the transition 
is similar over the different protocols. 



B. Isostatic and geometrical coordination numbers and symmetry breaking 



While the isostatic coordination number has been well documented at RCP [4],ll4|,ll8], the 
possibility of states with Zj = 6 along the coexistence region with rcp < <pj < me it requires 
more elaboration [6| . We recall, however that the limiting condition Z lso = 6 is necessary, but 
not sufficient, for a rigid isostatic aggregate: a kinematic condition for rigidity must hold 
[3] where the functions that describe the connections between the centers of contacting 
particles are independent. The presence of crystal-like regions suggests that this condition 
may not be satisfied. Thus, the packings with Zj = 6 in the coexistence region are not 
necessarily isostatic, except exactly at the freezing point. 

It is interesting to understand the geometrical rearrangements occurring during the RCP 
transition in light of the fact that the packings enter the coexistence region from RCP by 
keeping Zj = 6 constant. The particles modify the positions of the 6 contacts in average at 
RCP to create crystal-like regions without creating new contacts or destroying old ones. This 
implies that the arrangements of particles are such that particles in the second coordination 
shell come closer to the central particle and contribute to the first coordination shell, yet 
without producing a new contacting force since Zj is kept at 6 in the entire coexistence 
region. The new particles moving into the first coordination shell can be considered in 
geometrical contact but carrying no force. Thus following [7] we introduce the idea of 
geometrical contact, z g , as those particles in the first coordination shell that do not provide 
any force but still are close enough to the central particle to contribute to the geometrical 
contact network. The geometrical coordination number z g is different from the mechanical 
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coordination number Zj which only counts those contacts with nonzero forces. By definition 

Zg > Zj. 

While Zj is easy to measure in computer simulations of soft particles as the number of 
contacts between overlapping particles, the geometrical coordination number, z g , can be 
measured by slightly inflating the spheres up to 4% of their diameters and counting the 
resulting contact particles, as discussed in [7|. In practice, the geometrical coordination 
number measures the particles surrounding a central one with a gap between them from 
zero or negative (giving Zj) up to 5 = 0.08R as discussed in Q]. We notice that 5 = 0.08-R is 
much smaller than the location of the second peak in the radial distribution function which 
occurs around 5 ~ 2R. The value 5 = 0.08i? is specific for a system of iV=10,000. We 
expect this value to diminish in the thermodynamic limit. 

Figure 03 plots the geometrical and mechanical coordination z g and Zj as a function of 
(bj for the same packings as in Fig. HJ We find that along the disordered branch, z„ ~ Z~ 

n 

as expected [7|. However, in the coexistence region, Zj = 6 stays constant while z g keeps 
growing with <ftj. The separation between z g and Zj is a signature of the onset of ordering at 
the freezing point. As explained above, the system starts to crystallize by allowing particles 
in the second coordination shell to come closer to the central particle and moving the Zj 
contacting particles towards a FCC arrangement. At the melting point, the condition Zj = 6 
cannot hold any longer and the system transitions to the other branch with an increase of 
Zj up to 12. 

The distinction between z g and Zj is not only important for a characterization of the 
transition. It is also important to interpret the experimental results. Due to the uncertainty 
in detecting the exact position of the particles in any experiment, the exact mechanical 
coordination might be very difficult to obtain. Thus, a small uncertainty in the determination 
of the contacting particles 5 = 0.08.R will produce z g as shown in Fig. [5j3. One way to obtain 
the actual mechanical coordination from experimental data is to use the experimentally 
obtained coordinates of the particles as initial positions of a MD simulation using Hertz- 
Mindlin forces to relax the configurations and find the exact force balance network of the 
packing. Codes to develop this procedure are available at www.jamlab.org. We also provide 
most of the packings used in this study as well as the code to calculate the entropy. 
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C. Relation to the glass transition 



The thermodynamic character of the RCP transition seemingly contrasts to the non- 
thermodynamic viewpoint of the glass transition which proposes a dynamical arrest upon 
supercooling [l^]. However, the same phenomenology of vitrification could be applied to 
the RCP transition by extrapolating the entropy of the disordered branch s(4>j) into a 
metastable region below the freezing point as schematically shown in Figs. [1] and 03 from 
a — > b. Two scenarios may emerge: the metastable branch may end in a metastability limit 
at the spinodal dX/dS = [17( or it may continue until the entropy of the metastable liquid 
is zero as shown in Fig. |3J3. Such a scenario would predict a Kauzmann density 4>k at which 
the entropy of the disordered branch vanishes, signaling the existence of an ideal jammed 
glass analogous to the Kauzmann temperature in glasses 

III. DESCRIPTIVE VIEWPOINT OF THE RCP TRANSITION 
A. Orientational order parameter 

The orientational order is measured by associating a set of spherical harmonics with every 



bond joining a sphere and its neighbors 19j: 



Qim(f) = Y lm (e(r), </>{?)), (6) 

where the {Yi m (9,4>)} are spherical harmonics, and 9{r) and <fi(r) are the polar angles of 
the bond. The local orientational order parameter qi for a particle % is given by rotationally 
invariant combinations of Qi m , 



1/2 

47r -Vlo I 2 



21 , 

m=—l 



(7) 



where Q lm i is averaged over N t neighbors of this particle, 



^m,i = ]y X>m(ry. (8) 
1 3=1 

We also consider the global orientational order Qi, as well as the third-order invariants 
Wj, which are 
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where the average is taken over all the Nf, bonds in the packing, 



Qlr 



W Qhn (r) , 



and 



W, 



E 



where the coefficients 



rfii,m,2,m3 
mi+m 2 +m,3=0 
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bonds 
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(10) 



(11) 



are the Wigner 3j symbols. 

mi m>2 ^3 ' 

We use a definition of bond as in [19[ where all spheres within r c = 1.2c? of a given sphere 
are near neighbors, where d = 2R is the sphere diameter. We note that r c = 1.2d is the 
center of the first peak (centered at r = 0) and the second peak (centered at r = 1.4c?) 
in the radial distribution function g(r) of a perfect FCC lattice. Thus, this criteria is also 
consistent with the definition used in T. M. Truskett, S. Torquato, P. G. Debenedetti, Phys. 
Rev. E 62, 993 (2000), where r c is the first minimum in g(r). The neighbors can also 
be defined as those who have mechanical contacts with the given sphere, Zj. The basic 
results do not change by changing the definition of nearest neighbors. The results of the 
global orientational order Qi and Wj are shown in Fig. |2JA. and Fig. [2j3. To investigate the 
structures of the phases, Fig. |2C plots the probability distribution P(qe), which is the most 
sensitive measure of the local order parameters. 



B. Orientational correlation function 



The bond-angle correlation functions (Fig. [6]) can be obtained via 19j : 



Gi (r) = — (Qlrn (r) Q lm (0)}, (12) 

m=— I 

where the angular bracket indicates an average over all particles separated by f. A non-zero 
asymptotic value of Gq (r) implies a long-range correlation in the orientational order. From 
Fig. [6l it is clear that the crystal phase has long-range orientational order, while no order 
can be found in the random phase. 
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FIG. 6: The orientational correlation function Eq. (| 12|) for different packings with cfij. G§{r) 
vanishes at large r below <pj = 0.64, while it approaches a nonzero constant when cpj > 0.64. The 
nonzero asymptote of G§{r) is a signature of long-range correlation of orientational order. RCP is 
a well defined singularity at 4>j = 0.64 where the orientational symmetry breaking occurs. 

C. Local orientational disorder 



We have established that crystalline structures appear in the coexistence region and 
solid-like branch. The remaining question is what kind of lattice structure dominates in the 
crystallized packings. Since the differences of the orientational order parameters, especially 
Qq, are not significant between different perfect lattice clusters (such as the icosahedral, 



FCC, HCP, BCC, and SC clusters, see 



19], Fig. 2) we apply another measure, the local 



orientational disorder, to identify the crystalline nuclei as defined in M. Bargiel and E. M. 
Tory, Adv. Powder Technol. 12, 533 (2001). 

For a given sphere i, let 9^\. be the angle between the jth and fcth neighbors. Furthermore 
let 9^, 6^ p , 9]jj° B be similarly calculated angles for the perfect 13-sphere fragments of FCC, 



HCP and icosahedral packings, 9^ c be the angles for the 8-sphere fragment of a perfect BCC 
packing, and 9 s f k be the angles for the 6-sphere fragment of a perfect SC packing. The local 
disorders are defined as following: 



T fee 



11 12 



j=i k=j+i 



11 12 



\J 3=1 k=j+l 



(13) 



(14) 
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FIG. 7: (A) Distributions of local disorder #J CC denned in Eq. (fT3l) for packings with different 
(j>j as indicated. The distribution functions below <pj = 0.64 are Gaussian, the center of the peak 
moves to 6f cc = 0.2 as the volume fraction <j)j approaches 0.64. This disordered peak decreases 
above <j>j = 0.64 and seems to disappear above 4>j = 0.68. Two ordered peaks appear after RCP, 
the one at #{ cc = corresponds to FCC structure, while the other one at 0j €C = 0.16 corresponds to 
HCP. The peak at zero eventually evolves to a delta function as the packing structure approaches 
a perfect FCC. We use the median of the FCC peak and the disordered peak at <pj = 0.64 as a 
cutoff to identify local FCC structure, ie., particles with #J CC < 6^ c are defined as FCC crystalline 
nuclei, where 6^ c = 0.1. (B) Distributions of local disorder 9f cp defined in Eq. (j!4[) for packings 
with different as indicated. The distribution functions have similar behavior as those of #J CC . 
The cutoff 6c Cp = 0.075 is used to identify HCP crystalline nuclei. 
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Note that to calculate the values properly, we first need to sort the angles 9ij k , 9^, 6*^ p , 



#}fc OS > #jfc C an d Q B jk-> an d then compare them one by one. Also note that the BCC and SC 



3icos nhcc 
jk J °jk 

clusters have fewer neighbors than FCC, HCP and icosahedral clusters. 
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Since the 0j's measure the local disorder in the packings compared to a particular lattice 
structure, a perfect lattice cluster would have a zero value of 9. Any packing with significant 
amount of certain lattice clusters would indicate a peak centered around the origin in the 
distribution function of the local disorder 9i corresponding to that particular lattice. 

The distribution of FCC clusters P(9f c ) for packings with different <pj is shown in Fig. 
[7JA.. We find that FCC and HCP dominate in the crystalline packings for (pj > </> me it- Indeed 
we observe two prominent peaks in the distribution, one at FCC 9f c = and the other at 
HCP 9f c = 0.16, while BCC, SC and icosahedral ordering are negligible. The FCC peak 
dominance indicates that the majority of the clusters are FCC with a small proportion 
of HCP clusters. Similar to the distributions of local orientational orders shown in Fig. 
[2C, we find no crystalline clusters in the random packings as evidenced by the Gaussian 
distributions of P(9 i j cc ) for <pj < <ft Tcp as seen in Fig. [7JA.. In the coexistence region, the 
distributions are formed by a linear combination of different phases at melting and freezing. 
The distribution of HCP clusters P(# 4 hcp ) has similar behavior as P{9f c ) ) as showed in Fig. 



D. Cluster analysis of crystalline regions and correlation length 

We are able to define crystalline or nearly crystalline clusters in a packing based on the 
local orientational disorder and mechanical contacts, and visualize them in a 3d plot. The 
clusters are defined as follows: First, each node in the clusters is a sphere with 9f c < 9^ c 
or 9^ cp < 9^ cp , where 6^ cc = 0.1 and 6£ cp = 0.075, as determined in Fig. [7J The definition 
ensures that the first peak in P(9f c ) at zero in the distribution function of Fig. [7J\. is 
included in this consideration (as well as the analogous analysis for HCP). Next, if any two 
nodes are in mechanical contact, we build a link between the two nodes. Then the crystalline 
clusters are those nodes that are linked together. The clusters are visualized in Fig. [TJ 

Based on the definition of the crystalline clusters, the size of the largest cluster in the 
system and the correlation length of the clusters, £, are measured near the melting point. 
To calculate the correlation length, we first introduce the radius of gyration, R g (s), of a 
cluster consisting of s particles: 



EB. 




(18) 
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FIG. 8: Correlation length £ of crystalline clusters. To calculate the correlation length, we first 
introduce the radius of gyration, R g (s), of a cluster consisting of s particles. The correlation length 
of a perfect FCC lattice with periodic boundary condition is L/2, where L is the system size, so 
the correlation length is scaled by L/2 in the figure. The crystalline clusters start to percolate at 
4>j = 0.68 as £/(L/2) reaches a plateau with value 1. The results also show that the percolation at 
0meit does not depend on the system size. 

then the correlation length is given by 

2 s TR 2 Js)s 2 n< 



s 



^2 s2ris 



(19) 



where n s is the number of clusters of size s in the packing. 

The correlation length £ of the crystalline clusters is measured near the melting point. 
Figure [8] confirms the linear increase of the size of crystals along the coexistence region from 
the freezing point where £ = to the melting point. When the system melts at me it, £ 
reaches a plateau consistent with the system size. 



E. Radial distribution function 

The radial distribution function g(r) of the packing with volume fraction 0.72 in Fig. [9] 
shows all the peaks in FCC and HCP packings, which are indications of long range spatial 
order. On the other hand, the radial distribution functions of random packings only have 
two or three peaks (the second peak splits at <f>j = 0.64), corresponding to short range order. 
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FIG. 9: Radial distribution functions g(r) with different volume fractions <j>j. The figure clearly 
show the increasing of long range spatial order above <j>j = 0.64. 

In the coexistence region, g(r) has more peaks than those of random packings, but the 
magnitude of the peaks decays very fast as the distance r becomes larger. The long-range 
order increases with the volume fraction in the coexistence region. 

IV. THERMODYNAMIC VIEWPOINT OF THE RCP TRANSITION 

A. Free energy 



In the thermodynamics of jammed matter, the internal energy U is replaced by the 
volume W (usually called the volume function) 2l| . Other thermodynamic potentials, such 
as enthalpy H, free energy F, Gibbs free energy G, are related to the volume function W as 



H = W } 
F = W-XS, 
G = F = W - XS, 



(20) 
(21) 
(22) 
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where X is the compactivity of the system and S the entropy, related as: 



- " — (23) 

Since thermodynamic pressure is not considered in our case, the volume W and enthalpy 
H are equivalent and the Helmholtz free energy F and Gibbs free energy G are identical, 
so that we refer to them only as the free energy. The differentials of the thermodynamic 
potentials are 



dW = XdS + fidN, (24) 
dF = -SdX + fidN, (25) 

where \i is the chemical potential and iV is the number of particles. 

The free energy F, or Gibbs free energy G is equal to the product of N and //, 

F = C = N/2. (26) 



The free energy and chemical potential are continuous in the first order phase transition. 

dfi 



On the other hand, the derivative of the chemical potential, S = — t£, is discontinuous in a 



first order phase transition. 

The compactivity X, as well as the entropy density s = S/N, can be calculated from 
the fluctuations of the Voronoi volumes, which is an analogy of the Edwards theory to 
the standard Boltzmann statistical mechanics. The definition of a Voronoi cell is a convex 
polygon whose interior consists of all points closer to a given particle than to any other. 

n 

The Voronoi volume of a given particle, V mr , is the volume of such a Voronoi cell (see |7[ for 
more details). 

We define the Voronoi fluctuations as a 2 = (w^ ox ) — (w vor ) 2 , where w vor = V vor /V g is the 
reduced Voronoi volume of each particle and the average is done over all the particles in a 
packing. The quantities <J\{<j)j) and o- 2 ((j)j) denote the Voronoi fluctuations in the disordered 
and ordered phases, respectively, function of 6-i- 

We use the Einstein fluctuation theorem for jammed matter which is obtained from the 
similar relation in equilibrium systems by replacing the energy by the volume fluctuations 



23j-|25j: 
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{{5Wf) = k B X^. (27) 

Here we assume that ks plays the role of the Boltzmann constant in thermodynamics. Its 
value could be set to unity without changing the obtained results since in this context it just 
defines the units of entropy (in the main text we set ks = 1 to simplify). This means that 
we measure the compactivity in units of volume V g and that the entropy, which has units 
of ks, is dimensionless. In terms of the Voronoi fluctuation a and volume fraction <ftj, Eq. 
(1271) reads: 

Equation (I28I) applies to the disordered and ordered branches, o\ and cr 2 , separately. It 
does not apply to the coexistence region, since it requires the pure phases to calculate the 
fluctuations. 

Using Eq. ( 1281) we can calculate the thermodynamic quantities by integration. We first 
obtain the compactivity from: 



k 
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v g J^ p x rlp 

1 k B (** dd) 1 



+ j 0rlp < <f>j < 0rcp, (29a) 



v . . - , ■ I , 2 2fJ t \+Y ' 0melt < 0i < 0fcc, (29b) 

where X melt = X(0 me i t ) is the compactivity of the packing at the melting point and X rlp = 
X(0 rlp ) at RLP. 

The entropy density is then obtained by a second integration using the definition Eq. 
f[23l) which reads: 



1 <\>) ds 



X V g dfa 



(30) 



We integrate Eq. (130]) for each branch to obtain: 



1 c P 

*('">;) = «i,- P + V :l 1^ X (^^ 2 ' n ' i!> ~ " j ~ " nv ' (31a) 

^fcc 

= V 9 I x(^ 2 ' ^ melt ~ ^ ~ ^ fcc ' (31i>) 
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The entropy of FCC, Sf cc , is zero in the thermodynamic limit. 

Equations ([29]) and (1311) require three constants of integration: X r \ p , X me \ t and the entropy 
of RCP: s rcp . We now introduce three extra constraints to close the system. First, there are 
two conditions for equilibrium between two phases in jammed matter 17|, |2l|: (a) "thermal" 
equilibrium 

-Xmeit = X Icp = X c , (32) 

where X c is the critical compactivity at the transition, (b) The equality of the chemical 
potentials of the two phases at the melting and the freezing RCP point: /i me it = A^reez due to 
the conservation of number of particles. This is analogous to the equality of the free energy 
density at the transition from Eq. (T2T)]) . which allows to calculate the entropy at the freezing 
point via / me it = f rcp - 

s rcp — s melt T Vg y J, \3<J) 



and the entropy of fusion is then: 

As fus = S rcp - S me it = V g ^ rc P ^melt ^ _ ^ 

The final constant of integration to be obtained is the compactivity at the RLP point. 
As a first order approximation, X r i p can^bej.aken as infinite as has been shown in previous 



experimental and numerical studies 23|, l25| . However, when we compare the obtained 



entropy Eq. ( l3Ti) with an independent measure of the entropy using Shannon information 
theory (explained below) we find that there are slightly differences between both values of 
the entropy. Therefore, we consider X v \ p as a fitting parameter to be obtained by fitting 
the result of Eq. fl3T|) with the Shannon entropy which in principle does not require any 
integration constant to be calculated. We note that using the fitted value of X rlp instead 
of infinity does not change the final results, specifically the values of X c and the entropy of 
fusion, even though the obtained X T \ P is "far" from infinite. 

We summarize the calculation as follows: (a) We assume a value of X r \ p (which is later 
fitted with Shannon entropy) and integrate Eq. (I29al) from RLP to <p Tcp and obtain the 
compactivity X rcp = X(0 rcp ). (b) Using Eq. (132]) we obtain X me i t (or X c ). (c) Using X me i t , 
we integrate Eq. f]29bj) to obtain the X{<pj) in the ordered branch, thus completing the 
calculation of the compactivity equation of state, (d) We integrate Eq. (I31bl) up to mc it to 
obtain the entropy of the melting point: s me it- (e) Using Eq. ( |33l) we obtain the entropy at 
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RCP, s rcp . (f) This value is then plugged into Eq. (13 lap to finish the calculation of s(4>j) 
by a final integration. The above procedure is repeated for different values of X r i p starting 
from infinite up to a finite value that will match the Shannon entropy as discussed below. 

B. Shannon entropy 

Entropy of jammed matter can be calculated in two different ways: (i) The entropy 
from fluctuation theory explained above, (ii) The entropy from information theory [^J, so 
called " Shannon entropy" , related to configurational disorder since topologically equivalent 
structures are considered as the same state. Shannon entropy attempts to measure the 
disorder in a string of information as defined in the seminal work of Shannon. This concept 
has been adapted to the measurement of the configurational entropy in physical systems 
defined through a contact network by Vink and Barkema [R. L. C. Vink, G. T. Barkema, 
Phys. Rev. Lett. 89, 076405 (2002)]. In this work it was shown that the entropy obtained 
from the thermodynamic integration of fluctuations (or equivalently the specific heat) and 
the Shannon entropy are equivalent and accurately describe amorphous silicon and vitreous 
silica networks. The method has been extended to calculate the entropy of packings of 



granular materials in [23J. Below we explain the main details. 

The advantage of the Shannon entropy calculation over the thermodynamic integration 
is that it does not require a constant of integration as in Eq. ( 13~TT) . For each state of jammed 
matter, we associate a probability of occurrence Pi to the state, which is calculated as follows. 

We use the Voronoi cell and Delaunay triangulation for each particle to define a Voronoi 
network by considering contacts when a Voronoi side is shared between two particles, and 
hence are Delaunay contacts. A graph is constructed as a cluster of n particles that are 
Delaunay contacts, and by means of graph automorphism [B. D. McKay, Nauty user's guide 
(version 1.5), Tech. Rep. TR-CS-90-02, Australian National University (1990)] can be 
transformed into a standard form or "class" i of topologically equivalent graphs with a 
probability of occurrence p(i). In practice, we determine p(i) by extracting a large number 
m of clusters of size n from the system and count the number of times, /», a cluster i is 
observed, such that: 

p(i) = fi/m. (35) 
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Then the Shannon entropy is defined as: 

H{n) = -y^pjlapi, (36) 

where we have again assumed the "Boltzmann-like" constant in front to be one. The Shannon 
entropy density is obtained as: 

s shan = lim [H{n + 1)-H {n)\ , (37) 

by linear fitting of the extensive part of the Shannon entropy 

In general, the fluctuation entropy from Eq. (|3TT) is greater than the Shannon entropy 
from Eq. (I3T1) because Shannon entropy only counts configurational disorder and additional 
entropy could arise from freedom to move grains within the clusters of n particles without 
disrupting the Delaunay network. However, we discover that the fluctuation entropy ob- 
tained from Eq. fl3Tl) and Shannon entropy from Eq. (1371) only differ by a scaling constant 
k = 0.1, ie., ks = s s h an , see Fig. [3j3. Beyond this multiplicative constant the agreement 
between both estimations of the entropy is very good. Due to finite size effects the Shannon 
entropy also gives a nonzero value of the entropy of FCC, Sf cc = 0.6. This value is subtracted 
from the calculations. By fitting s s h an with the thermodynamic entropy we obtain the final 
constant of integration in Eq. (I29al) . X r i p = 0.51^. While it is obvious that this constant 



is far from the infinite value which is expected and used in [23j, |25| for the RLP limit, we 
notice that our final results are not very sensitive to the exact value of X r i p . For instance, 
by setting X v \ p — > oo the fitting of the fluctuation entropy is slightly off in comparison with 
Sshan only in the vicinity of RLP but the values of X c and the entropy of fusion do not have 
appreciable change. 

Figure H] in the main text shows the thermodynamic quantities in the first order phase 
transition of jammed matter. They are consistent with the general thermodynamic picture. 
Figure [10] plots the entropy versus w showing a linear dependence in the coexistence region. 
The entropy is an interpolation of the form: s x = xs m eit + (1 — x)sf reez where x is the 
concentration of crystal clusters in the coexistence. Since X = dW/dS, the linearity of s 
between 0.64 and 0.68 is a manifestation of the coexistence of two phases at a constant X c . 
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